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1.  SUMMARY 


1.1.  Objective  and  Accomplishments 

The  overall  objective  of  the  proposed  research  was  to  provide  a  quantitative  description  of 
the  microstructurd,  as  well  as  macroscopic  mechanical  behavior  of  particulate  materials  with 
viscous  and  viscoelastic  intergranular  cement.  By  cement,  a  material  is  meant  which  fills  the 
space  between  two  or  more  particles  at  their  contacts.  The  volumetric  content  of  cement  in 
the  pore  space  of  an  aggregate  may  vary  from  zero  to  100%. 

This  objective  was  achieved  by  deriving  microstructural  contact  laws  for  the  combination 
of  two  elastic  spherical  grains  with  viscous  cement.  A  simple  and  accurate  viscoelastic 
approximation  (Cole-Cole  model)  was  found  for  the  rigorous  micromechanical  solution.  The 
contact  law  thus  obtained  was  used  to  calculate  the  effective  elastic  moduli  of  a  granular 
aggregate  with  viscous  cement  at  a  varying  frequency  of  oscillations.  The  theoretical  results 
were  verified  by  experiments  conducted  on  a  pack  of  glass  beads  with  viscous  epoxy  cement. 

Our  second  achievement  was  the  extension  of  an  exact  solution  for  the  elastic  moduli  of  a 
cemented  grain  pack  at  a  low  cement  content  to  a  solution  for  the  case  where  cement  occupies 
the  entire  pore  space  or  a  large  portion  of  it.  The  theoretical  estimates  accurately  mimic  the 
experimental  elastic-wave  velocity  measurements  on  artificial  and  natural  granular  cemented 
materials. 

The  main  relevance  to  the  Air  Force  mission  is  through  a  quantitative  description  of 
cemented  geomaterials  such  as  asphalt  cement.  The  results  have  been  used  at  the  Wright 
Laboratory  to  model  asphalt  concrete's  behavior. 

1.2.  Personnel  Supported 

Senior  Personnel:  Dvorkin,  Jack,  PI  (US).  Ph.D.  Students:  Ran  Bachrach  (non-US);  Amy 
Day-Lewis  (US);  Doron  Galmudi  (non-US);  Mike  Helgerud  (US);  Yuguang  Liu  (non-US); 
James  Packwood  (US);  Li  Teng  (non-US). 

1.3.  Publications 

1.  Jack  Dvorkin,  Mickaele  Le  Ravalec,  Jim  Berryman,  and  Amos  Nur,  1997,  Effective 
moduli  of  particulates  with  elastic  cement,  submitted  to  Mechanics  of  Materials. 

2.  Klaus  C.  Leurer  and  Jack  Dvorkin,  1997,  Intergranular  Squirt  Flow  in  Sand:  Grains 
with  Viseous  Cement,  submitted  to  Int  J.  Solids.  Struct. 

1.4.  Interactions  and  Transitions 
Workshops  and  Conferences: 

•  1st  Stanford  Annual  Gas  Hydrate  Workshop;  Stanford,  CA,  1997:  Modeling  Hvdrate- 
Cemented  Sands. 

•  1997  Stanford  Rock  Physics  Meeting;  Stanford,  CA,  1997:  Natural  Cemented  Rocks. 

•  Invited  Lecture  at  Mobil  Technical  center  in  Dallas,  1997:  Cemented  Rocks. 

•  McNU  1997  Conference:  Transition  From  High-Porositv  to  Filled  Particulates:  Effective 
Medium  Modeling. 

Transitions: 

•  USAF  Laboratories:  Wright  Laboratory,  Tyndall  AFB,  Dr.  Han  Zhu  —  used  the 
cementation  theory  (viscoelastic  cement)  to  describe  the  deformation  of  asphalt  concrete. 

1.6.  Patents  AND  Inventions 
No  patents  or  inventions  resulted  form  this  effort. 
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2.  TECHNICAL  PART 


In  this  part  the  technical  results  are  presented  in  the  form  of  journal  articles.  The 
numeration  of  formulas  and  figures  is  internal  for  every  section. 

2.1.  Intergranular  Squirt  Flow  in  Sand:  Viscous  Cement 

SUPPORTED  IN  PART  BY  THE  NAVAL  RESEARCH  LABORATORY 

ABSTRACT 

We  obtained  an  exact  solution  to  the  problem  of  deformation  of  two  elastic  spherical 
particles  with  viscous  cement  at  their  contact.  This  model  is  intended  to  mimic  granular 
marine  sediment  whose  grains  are  covered  with  viscous  fluid  of  geologic  (clay  suspension) 
or  biogenic  nature.  It  is  generally  applicable  to  a  granular  aggregate  with  viscous  cement. 
The  solution  for  the  normal  stiffness  of  a  two-grain  combination  is  reduced  to  an  ordinary 
integro-differential  equation  that  has  to  be  solved  numerically.  Based  on  the  numerical 
solution,  we  found  two  approximate  expressions  for  the  normal  stiffness  of  a  two-grain 
combination.  The  first  approximation  based  on  the  Maxwell  viscoelastic  model  is  not  very 
accurate.  The  second  approximation  based  on  the  Cole-Cole  model  is  very  accurate  and 
can  be  used  instead  of  the  numerical  solution.  The  effective  elastic  moduli  of  the  aggregate 
can  be  calculated  using  statistical  averaging  for  a  dense  random  pack  of  identical  spheres. 
The  theoretical  results  match  well  experimental  data  obtained  on  a  glass  bead  pack  with 
viscous  epoxy  cement. 

INTRODUCTION  AND  PROBLEM  FORMULATION 

The  elastic  properties  of  a  granular  aggregate  such  as  oil  sand  or  marine  sediment 
strongly  depend  on  the  stiffness  of  the  grain-to-grain  contacts.  The  stiffness  of  the  contact 
region  can  be  affected  by  the  presence  of  high- viscosity  fluid  (e.g.,  heavy  oil,  clay 
suspension,  or  a  biogenic  material)  that  envelops  the  grains.  This  fluid  may  act  as  contact 
cement  by  reinforcing  the  intergranular  contacts  (Figure  1). 
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Figure  1.  Left:  two  spherical  grains  enveloped  by  high- viscosity  fluid.  Right:  the  grain  contact 
region  with  viscous  contact  cement. 

The  effective  elastic  properties  of  a  granular  aggregate  with  viscous  cement  are 
frequency-dependent.  At  a  low  frequency,  hydrodynamic  pressure  in  viscous  cement  will 
be  at  equilibrium  with  the  surrounding  pore  space.  Therefore,  viscous  cement  will  not 
contribute  to  the  stiffness  of  the  grain-to-grain  contact.  At  a  very  high  frequency,  viscous 
cement  becomes  unrelaxed  and  deforms  as  an  elastic  body.  Now  it  may  strongly  reinforce 
the  contacts.  At  intermediate  frequencies  local  (squirt)  flow  of  viscous  cement  develops 
between  the  grains.  This  viscous  flow  is  responsible  for  velocity-frequency  dispersion  and 
attenuation.  Our  goal  is  to  quantitatively  describe  the  squirt  flow  of  viscous  cement  and 
understand  how  it  affects  the  stiffness  of  the  sediment  frame  and,  therefore,  acoustic 
velocity  and  attenuation  in  the  sediment. 

In  order  to  solve  the  problem,  we  assume  that  the  grains  are  spherical  and  elastic,  and 
cement  is  a  Newtonian  compressible  viscous  fluid.  The  main  part  of  the  solution  is  finding 
the  normal  stiffness  of  a  two-grain  combination.  This  normal  stiffness  is  defined  as  the 
ratio  of  the  applied  (to  the  grains)  normal  force  to  the  resulting  increment  of  the 
displacement  of  the  sphere  center.  We  assume  that  the  shear  stiffness  of  this  combination 
is  non-existent.  Also,  we  assume  that  the  effective  pressure  acting  on  the  grains  is  zero, 
i.e.,  the  direct  grain-to-grain  contact  is  a  point  and  the  grains  are  in  close-pack  suspension. 
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Once  the  expression  for  the  normal  contact  stiffness  is  obtained,  it  can  be  used  in  a 
numerical  discrete  element  code.  For  a  special  case  of  a  random  close  pack  of  identical 
spherical  grains  the  following  analytical  expressions  can  be  used  to  calculate  the  effective 
bulk  and  shear  (G^^)  moduli  of  the  aggregate  (e.g.,  Dvorkin  et  al,  1994): 


^  _n(l-0)p  ^  . 

\1kR 


^Eff 


(1) 


where  is  the  normal  stiffness  between  two  grains;  n  ~  9  is  the  coordination  number  (the 
average  number  of  contacts  per  grain);  0  «  0.36  is  the  aggregate's  porosity. 


GOVERNING  EQUATION 

We  describe  the  dynamics  of  the  viscous  cement  by  examining  its  viscous  flow  induced 
by  the  oscillations  of  the  grain  surfaces  (Dvorkin  et  al.,  1990).  Specifically,  we  consider 
oscillations  at  a  fixed  angular  frequency  (O.  Then  the  thickness  b  of  the  cement  layer  in 
the  vicinity  of  the  point  grain-to-grain  contact  is  a  function  of  the  radial  coordinate  r  and 
time  t: 

Kr,t)  =  -^  +  K(r)e‘",  (2) 

where  R  is  the  grain  radius,  and  «  b  is  the  amplitude  of  the  grain  surface  oscillations. 

For  an  acoustic  compressible  fluid  with  the  speed  of  sound  Cq  the  pressure  increment 
dP  is  proportional  to  the  density  increment  dp  as 

dP  =  CQ^dp.  (3) 

The  mass  conservation  equation  in  the  contact  region  (Figure  1)  is 

dp  ^  d{pu)  ^  pu  ^  d(pw)  ^ 
dt  dr  r  dz 
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where  u  is  the  velocity  component  along  the  r  coordinate,  and  w  is  the  velocity 
component  along  the  z  coordinate. 

Let  us  integrate  equation  (4)  in  the  z  direction  from  0  to  and  take  into  account  that 


dJy 

w  =  0,  z  =  0;  w  =  —,z  =  b. 

at 


(5) 


Then 


d{bp)  ,  \,d{pu)  ,  pu 
dt 


J  dr  ^ 


(6) 


The  approximate  Navier-Stokes  equation  for  the  radial  flow  in  the  contact  gap  is 

du  1  dP  d^u 

dt  p  or  dz 


where  77  is  the  kinematic  viscosity.  A  solution  of  equation  (7)  in  the  frequency  domain  is 

(8) 


^  \  dP ..  cosh(z^jft)/  77) 

u{r,z)  =  — : - ^[1 - , /.  ■  ']• 

mp  dr  cosh(oym)/  77) 


We  assume  that  viscous  cement  is  acoustic  fluid  and  thus  variations  of  density  are  small 
as  compared  to  its  reference  value.  Then,  and  using  equations  (2)  and  (3),  we  have 


C/t  Cq  (yt 


dr 


dr  dt 


(9) 


Finally,  we  substitute  equation  (8)  into  (6)  and  use  equations  (2)  and  (9)  to  obtain  an 
ordinary  differential  equation  for  pressure  P  in  the  frequency  domain; 


( 


d^P  .  dP^,^  tanhA,  .  ^  dP  „  2  ^  ,..2  •  2/? 


dr^  rdr 


)(!■ 


)  +  2— -tanh  A  +  P— Y  = -ft)  pZ?o ■ 

rdr  Cn 


.2  ’ 


(10) 


X  =  (r^  /2R)^]iO)/  p. 

In  order  to  find  b^ir)  we  examine  the  elastic  deformation  of  the  grain  surface.  We  use 
the  following  compatibility  equation  among  b^ir),  the  displacement  V(r)  of  the  grain 
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surface  in  the  z  direction,  and  the  displacement  d  of  the  sphere  center  (Dvorkin  et  al, 
1994): 

h(r)  =  V(r)-6.  (11) 

Then  we  use  an  integral  equation  that  relates  y(r)  to  pressure  P{r)  in  the  viscous  flow  that 
is  exerted  upon  the  grain  surface  in  the  contact  region  (Timoshenko  and  Goodier,  1970): 


where  v  and  G  are  the  grain  Poisson's  ratio  and  shear  modulus,  respectively;  and  a  is  the 
radius  of  the  contact  region  (Figure  1).  This  radius  is  related  to  the  thickness  h  of  the  layer 
of  the  viscous  fluid  around  the  grain  as 

a  =  V^.  (13) 

The  integration  domain  in  equation  (12)  is  shown  in  Figure  2. 


Figure  2.  Integration  domain  in  equation  (12). 


Now  by  combining  equations  (10),  (1 1),  and  (12),  we  arrive  at  the  governing  equation 
for  pressure  P : 


The  boundary  conditions  for  equations  (10)  and  (14)  are  zero  pressure  fluctuation  (from  the 
ambient  hydrostatic  pressure)  at  the  open  boundary  of  the  viscous  cement  layer,  and  no 
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flow  at  the  center  of  the  layer.  The  expressions  for  these  boundary  conditions  are, 
respectively: 

P  =  0,r  =  a;  —  =  0,r  =  0.  (15) 

dr 


NUMERICAL  SOLUTION 


To  solve  equation  (14)  numerically  we  first  normalize  it  by  introducing  the  following 
notations: 


/  = 


P  f.  s  r  a  Cr.  .  iR^W  ^  2d 

-^,d  =  —,x  =  —,a  =  —,Y  =  —^,A= — . - ,  C  =  — . 

^  "  R  '  RO)  2  V  r?  R 


pc^  R  R 
The  resulting  equation  is: 


2,,  tanh/l,  2  2  /I  tanhA.  2r 

7  (1 - 7—)^  ttt  +  r  xa - ^  +  2 tanh"  A)^  +  x7  = 

A  ox  A  ox 


2/1  \  ^  ;tcos(0+'\/a^^^T^lin^|£ _ 

— P£oJ - —  jdcp  J  fisjx^  +  - 2jc<^  cos (p)d^  +  C. 

0  0 


(16) 


In  order  to  find  the  desired  normal  stiffness  between  two  grains  we  relate  the  force 
acting  on  the  grain  to  displacement  6.  Then 


a 


AttRpcq 


0 


(17) 


It  is  clear  now  that  equation  (16)  can  be  solved  with  an  arbitrarily  chosen  non-zero  C.  The 
resulting  normal  stiffness  will  not  depend  on  this  choice.  We  solve  equation  (16)  using  the 
quadrature  method  (e.g..  Delves  and  Mohamed,  1985). 

The  results  of  solving  equations  (16)  and  (17)  with  input  parameters  G  =  45  GPa;  v  = 
0.064;  Cq  =  1500  m/s;  p  =  1  g/cc;  R  =  10-4  m;  /j  =  io-7  m;  w  =  9;  ^  =  0.36;  and  cement 
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viscosity  1  cPs;  10^  cPs;  and  10^  cPs  are  given  (in  terms  of  the  real  and  imaginary  parts  of 
the  normal  stiffness  SJ  in  Figure  3.  The  transition  zone  (from  the  low-frequency  to  high- 
frequency  limit)  of  the  real  part  of  S„  as  well  as  the  peak  of  the  imaginary  part  of  move 
up  the  frequency  axis  as  viscosity  of  cement  decreases.  The  frequency  of  this  transition  is 
inversely  proportional  to  cement  viscosity. 


Figure  3.  Real  and  imaginary  parts  of  the  normal  stiffness  versus  frequency.  The  viscosity  of  cement 
is  given  in  the  plots. 

Once  S„  is  available,  we  can  calculate  the  effective  bulk  and  shear  moduli  of  the 
aggregate  using  equations  (1).  We  stress  that  equations  (1)  have  been  derived  for  the  case 
of  static  deformation  of  a  granular  aggregate.  In  the  case  under  consideration  the  response 
of  the  aggregate  to  the  load  is  clearly  time-dependent  (Figure  3).  Moreover,  the  imaginary 
part  of  this  response  has  the  same  order  of  magnitude  as  the  real  part.  This  fact  calls  for  a 
special  non-static  analysis  of  this  case.  We  will  address  this  issue  in  future  work.  As  of 
now,  we  assume  that  equations  (1)  are  applicable  to  the  case  under  examination. 
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IMPLICATIONS  OF  SOLUTION 


We  use  equations  (1)  to  find  the  effective  bulk  and  shear  moduli  of  the  cemented 
aggregate,  and  then  equations 


U  = 


'=1 


EffJ 


Pm 


Re(<j£^) 

i  Pm 


(18) 


where  is  the  bulk  density  of  the  aggregate,  to  calculate  compressional-  and  shear-wave 
velocities  and  V^.  The  grain  density  in  these  calculations  was  2.65  g/cc  (quartz). 


Figure  4.  P-  and  S-wave  velocity  versus  frequency  for  dry  frame  and  saturated  sediment. 


Equation  (18)  gives  wave  velocities  for  the  dry  cemented  frame  of  marine  sediment.  In 
order  to  calculate  these  velocities  in  the  water-saturated  sediment,  we  use  Gassmann’s 
(1951)  equation  whose  applicability  to  the  case  under  examination  (see  discussion  at  the 
end  of  the  previous  section)  will  be  investigated  later.  The  result  for  cement  viscosity  10“^ 
cPs  and  the  saturating  fluid  being  pure  water  are  given  in  Figure  4.  The  blow-ups  for  P- 
wave  velocity  are  given  in  Figure  5. 
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Figure  5.  P-  wave  velocity  versus  frequency  for  dry  frame  and  saturated  sediment  (blowup). 


APPROXIMATE  SOLUTION:  MAXWELL  BODY 

By  carrying  out  the  exact  numerical  solution  for  the  normal  stiffness  in  a  wide  range  of 
input  parameters,  we  find  the  following  approximate  expressions: 

Re(5J  ^  .  Im(5„)  _  . 

2kRpcq  “  1  +  ’  2nRpc^  “  1  +  ’ 

where 

=  Aa^  +  Ba  +  C\ 

A  =  19.7472  •  exp(-5 1.9238 A),  B  =  0.795673  •  A"® 

C  =  -0.00391064 •  A  =  pc^\\ -  v)  /  {nG)-,  x  =  0.8^27] / 


(19) 


(20) 


Equations  (19)  describe  a  Maxwell  body  (e.g.,  Bourbie  et  al.,  1987)  with  the  following 
constitutive  equation: 


dt  p  E  dt’ 


(21) 


where  F  is  the  contact  force  and  6  is  the  corresponding  displacement;  and  viscoelastic 
constants  p  and  E  are: 
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E  =  2KRpc^S^,  P  =  2kRpcqS„t.  (22) 

This  approximate  solution  is  compared  to  the  exact  numerical  solution  in  Figure  6.  The 


high-frequency  and  low-frequency  end  members  of  the  exact  and  approximate  solutions  are 
very  close  to  each  other,  however,  the  amplitude  of  the  imaginary  part  of  the  approximate 
solution  is  about  twice  that  of  the  exact  solution.  The  reason  is  that  in  the  Maxwell  body 
model  we  use  a  single  relaxation  time  which  results  in  a  narrow  (in  the  frequency  domain) 
transition  region  from  the  low-frequency  behavior  to  the  high-frequency  behavior. 
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Figure  6.  Real  and  imaginary  parts  of  the  normal  stiffness  versus  frequency,  exact  and  approximate 
Maxwell  solutions.  The  viscosity  of  cement  is  given  in  the  plots. 


APPROXIMATE  SOLUTION:  COLE-COLE  MODEL 


A  much  better  (than  the  Maxwell  body)  approximation  to  our  exact  numerical  solution 
is  given  by  the  Cole-Cole  (1941)  model: 


Re(5j 

2nRpcQ  “  1  +  •V2<ot  -h  m  ’ 

Im(5J  VfflT/2 

2kRpcq  “  \  +  -^2m  +  m' 


(23) 
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The  results  of  this  approximation  are  compared  with  the  exact  solution  in  Figure  7. 


02468  02468 

Log  Frequency  (Hz)  Log  Frequency  (Hz) 

Figure  7.  Real  and  imaginary  parts  of  the  normal  stiffness  versus  frequency,  exact  (thin  dark  curve)  and 
approximate  Cole-Cole  (bold  gray  curve)  solutions.  The  viscosity  of  cement  is  given  in  the  plots. 

EXPERIMENTAL  VERIFICATION  OF  RESULTS 

We  apply  our  theoretical  model  to  reproduce  ultrasonic  (1  MHz)  experimental  P-wave 
velocity  measurements  conducted  on  a  dense  random  pack  of  identical  glass  beads  partially 
saturated  with  liquid  epoxy.  The  measurements  have  been  conducted  at  finite  confining 
pressure  values  since  it  was  impossible  (due  to  coupling  problems)  to  propagate  a  pulse 
through  the  system  at  zero  confining  pressure.  Our  theoretical  model,  on  the  other  hand, 
assumes  that  the  grains  are  not  precompacted. 

To  apply  our  model  to  the  system  used  in  the  experiments,  we  modify  it  by  assuming 
that  an  elastic  spring  is  placed  in  parallel  to  the  original  viscoelastic  model.  This  spring 
represents  the  finite  stiffness  of  the  system  without  epoxy  present.  As  a  result,  we  can 
determine  the  stiffness  of  the  spring  using  velocity  measurements  at  zero  epoxy  saturation. 
Then,  in  order  to  convert  the  experimental  measurement  results  at  a  confining  pressure  to 
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those  without  confinement,  we  have  to  subtract  the  elastic  modulus  of  the  dry  system  from 
that  at  a  finite  epoxy  saturation. 


Figure  8.  Compressional-wave  modulus  (M -modulus)  difference  (saturated  minus  dry)  of  a  glass-bead 
pack  versus  epoxy  saturation.  Filled  circles  are  the  data,  triangles  are  from  our  model. 


We  have  available  velocity  measurements  5,  10, 15,  and  20  MPa  confining  pressure  in 
a  dry  glass  bead  pack  and  in  the  pack  with  25%  epoxy  saturation  (Yin,  1993).  In  order  to 
reproduce  the  experiment,  we  assume  that  the  shear  modulus,  Poisson's  ratio,  and  density 
of  glass  are  26.2  GPa,  0.277,  and  2.48  g/cm^,  respectively  (Dvorkin  et  al.,  1994). 
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We  did  not  have  direct  data  for  the  viscosity,  density,  and  the  acoustic  velocity  in  pure 
epoxy.  We  selected  reasonable  values  of  200  cPs,  1500  m/s,  and  1  g/cm^,  respectively 
(Yin,  1993).  These  values  allow  us  to  consistently  match  the  experimental  data  at  all 
confining  pressures  (Figure  8)  which  confirms  the  validity  of  our  theoretical  model. 

In  these  calculations,  we  assumed  that  the  epoxy  evenly  envelopes  every  grain. 
Therefore,  the  normalized  contact  radius  a  can  be  related  to  epoxy  saturation  as 
(Dvorkin  and  Nur,  1996) 


a  = 


2(1) 


3(1-0) 


(24) 


RELATING  ELASTIC  MODULI  TO  POROSITY 


The  above  results  are  appropriate  for  modeling  the  elastic  response  of  an 
unconsolidated  sediment  at  about  36%  porosity,  which  is  the  porosity  of  a  random  pack  of 
identical  spheres  (critical  porosity  ~  Nur  et  al.,  1991).  In  order  to  extend  this  result  for  the 
entire  porosity  range,  we  use  a  model  of  Dvorkin  et  al.  (1998)  where  the  Hashin-Shtrikman 
bounds  are  used  to  calculate  the  dry-frame  elastic  moduli  from  those  at  the  critical  porosity. 
For  porosity  <j)  below  critical  porosity  0^,  we  have 


^z,^(0)  =  [ 


0/0,  1-0/0.  .-1  4 

JL  AH-  1^  A.  A  fl  ^ 

^Ejf  ”  3  ^Ejf  TV  T  ^  J 


3  ^Ejf 


G^ff  +  Z  G  +  Z  6 


V  ^Bff  +  '^G^ff  j 


(25) 


0<0c’ 


where  and  come  from  equation  (1)  using  our  sphere  pack  model;  and  K  and  G 
are  the  bulk  and  shear  moduli  of  the  mineral  phase,  respectively.  For  porosity  above 
critical  porosity  we  have 
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(26) 


^Dry 

^Dry 


=  [ 
=  [ 


il-(^)lil-(t>,) 


l^Eff 

^  ((/>-f)/(l-^J 

z 


r'  __r 
J  2 

r'-z,  </>>-/»,. 


The  dry-frame  moduli  are  plotted  versus  porosity  in  Figure  9  where  the  parameters 
used  are  the  same  as  for  plots  in  Figure  4,  and  frequency  is  infinite  (elastic  limit). 


Figure  9.  Dry-frame  bulk  and  shear  moduli  versus  porosity  calculated  using  equation  (25)  and  (26). 

CONCLUSION 

The  theoretical  model  offered  here  is  a  first  step  towards  rigorously  describing  the 
visco-elastic  behavior  of  a  particulate  system  with  viscous  cement  at  grain  contacts.  Such 
systems  may  be  encountered  in  nature  (marine  sediments  with  clay  and/or  biogenic  viscous 
matter  enveloping  the  grains;  heavy  oil  sands;  shallow  sands  with  viscous  contaminant)  and 
in  engineering  applications  (asphalt  concrete).  The  model  developed  here  is  strictly 
speaking  appropriate  only  for  describing  suspensions  with  zero  effective  pressure.  Our 
next  step  is  to  extend  this  model  for  the  case  where  the  contacting  grains  are  subject  to  a 
finite  confining  force  and  thus  have  a  direct  contact  (Hertzian)  area  rather  than  a  point. 
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Such  theoretical  development  will  be  accompanied  by  pulse-transmission  experiments  to 
measure  elastic-wave  velocities  in  relevant  granular  systems. 

Another  goal  is  to  find  an  appropriate  physical  representation  for  the  Cole-Cole  model 
employed  here  that  can  be  easily  used  in  a  numerical  discrete  element  code. 

In  calculating  the  elastic  moduli  of  the  fully  saturated  system,  we  used  Gassmann's 
equation  where  the  dry-frame  moduli  are  given  by  our  model  for  a  pack  with  viscous 
cement.  We  did  so  for  simplicity.  However,  it  is  straightforward  to  use  the  Biot  model 
instead  which  may  help  account  for  the  global  frequency  dispersion  in  the  system. 
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2.2.  Effective  Moduli  of  Particulates  with  Elastic  Cement 


ABSTRACT 

We  give  an  approximate  method  for  calculating  the  effective  elastic  moduli  of  a  close 
random  pack  of  identical  elastic  spheres  whose  pore  space  is  partially  or  completely  filled 
with  elastic  cement.  To  construct  the  solution  we  start  with  the  uncemented  pack  and  add 
small  amounts  of  cement  uniformly  at  every  grain  contact.  The  effective  moduli  of  such  an 
aggregate  depend  on  those  of  the  grains  and  cement  (which  are,  generally,  different),  and 
on  the  amount  of  cement.  Moduli  for  grain/cement  mixtures  are  available  from  the  contact 
cementation  theory.  Next  we  introduce  an  isotropic  elastic  body  with  void  inclusions  (e.g., 
spherical)  whose  porosity  is  the  same  as  the  porosity  of  the  pack  with  contact  cement.  The 
moduli  of  the  matrix  are  as  yet  unknown.  We  find  them  by  assuming  that  the  effective 
moduli  of  the  elastic  body  with  inclusions  are  equal  to  those  of  the  pack  with  contact 
cement.  Then  an  appropriate  effective  medium  theory  (e.g.,  self-consistent  or  differential 
effective  medium  approximation)  provides  two  equations  for  the  bulk  and  shear  moduli  of 
the  matrix.  Finally,  we  use  the  chosen  effective  medium  theory  to  calculate  the  moduli  of 
the  elastic  body  with  the  matrix  thus  defined  whose  inclusions  are  now  filled  with  the 
cement.  We  assume  that  these  are  also  the  moduli  of  the  pack  whose  pore  space  is 
completely  filled  with  the  cement.  To  calculate  the  elastic  moduli  of  the  pack  whose  pore 
space  is  partially  filled  with  the  cement,  we  assume  that  they  are  those  of  the  elastic  body 
with  the  matrix  thus  defined  and  with  two  types  of  inclusions  —  void  and  filled  with 
cement.  The  porosities  of  this  body  and  of  the  cemented  pack  are  the  same.  This  result 
gives  one  solution  to  a  general  and  previously  unresolved  effective  medium  problem  of 
contacting  elastic  inclusions  embedded  in  an  elastic  matrix.  It  provides  a  transition  from 
high-porosity  cemented  to  completely  cement-filled  granular  materials.  The  solution 
accurately  predicts  the  compressional  elastic  modulus  (the  product  of  the  bulk  density  and 
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P-wave  velocity  squared)  from  wave-propagation  experiments  on  cemented  glass  beads 
and  natural  rocks. 

INTRODUCTION  AND  PROBLEM  FORMULATION 

When  estimating  the  effective  stiffness  of  many  composites,  concrete  and  natural  rock 
among  them,  one  often  faces  the  problem  of  calculating  the  moduli  of  an  elastic  matrix  with 
elastic  inclusions.  This  problem  has  been  extensively  addressed  in  the  past.  The  existing 
solutions  are  analyzed  and  summarized  by,  e.g.,  Christensen  (1991),  Zimmerman  (1991), 
Wang  and  Nur  (1992),  Nemat-Nasser  and  Hori  (1993),  and  Berryman  (1995).  Accurate 
solutions  are  available  for  small  inclusion  concentrations  and  various  inclusion  shapes.  To 
obtain  an  estimate  for  high  concentrations,  one  may  find  the  self-consistent  (SC)  or 
differential  effective  medium  (DEM)  approximations  to  be  useful  and  reasonably  accurate. 
Some  fundamental  solutions  for  composites  with  concentrated  interacting  inclusions  belong 
to  Frankel  and  Acrivos  (1967),  Goddard  (1977),  and  Chen  and  Acrivos  (1978). 

An  alternative  to  direct  estimates  is  the  bounding  of  the  overall  moduli  of  a  composite 
(Hashin  and  Shtrikman,  1963).  Nemat-Nasser  et  al.  (1993)  used  generalized  Hashin- 
Shtrikman  variational  principles  to  provide  accurate  bounds  for  composites  with  periodic 
microstructure.  They  also  introduced  a  "modified  equivalent  inclusion  method"  for  directly 
estimating  the  moduli  of  such  composites.  The  results  fall  between  the  upper  and  lower 
bounds  and  are  supported  by  experimental  data  on  the  shear  modulus  of  an  incompressible 
matrix  reinforced  by  rigid  spherical  particles. 

In  composites  with  concentrated  elastic  inclusions,  direct  interaction  between  the 
inclusions  through  the  thin  layers  of  the  matrix  is  very  important  (Goddard,  1977).  The 
contact  cementation  theory  (CCT)  of  Dvorkin  et  al.  (e.g.,  1994)  highlights  this  conclusion. 
It  allows  one  to  calculate  the  effective  elastic  properties  of  a  particulate  aggregate  where 
cement  is  located  around  the  grain  contacts  (Figure  la).  CCT  predicts  that  even  small 
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volumes  of  contact  cement  (which  may  be  soft  as  compared  to  the  grains)  act  to  strengthen 
the  aggregate's  elastic  moduli,  thus  turning  the  nonlinear  problem  of  contact  analysis  into  a 
simpler  linear  one  such  as  that  found  typically  in  composites  analysis.  The  initial  filling  of 
the  contact  gaps  is  most  important  —  the  effect  of  any  additional  cement  placed  around  this 
initial  filling  is  relatively  small.  This  result  is  supported  by  experiments  (Figure  2). 
Experiments  also  show  that  even  by  filling  the  entire  pore  space  of  a  sphere  pack  one 
cannot  achieve  the  drastic  relative  stiffness  increase  for  which  the  small  volumes  of  the 
contact  cement  are  responsible. 


Figure  1.  a.  Grains  with  contact  cement,  b.  Grains  with  pore  space  completely  filled,  c.  Contact 
cement  and  pore-filling  cement. 


The  assumptions  used  in  CCT  make  its  use  appropriate  only  for  high-porosity 
particulates  where  cement  is  located  in  the  immediate  vicinity  of  the  contacts.  As  is,  CCT 
cannot  be  used  to  estimate  the  elastic  constants  of  a  zero-  or  low-porosity  aggregate  where 
cement  fills  the  whole  pore  space  (Figure  lb)  or  large  portions  of  it. 

Therefore,  a  question  naturally  arises:  How  do  we  calculate  the  effective  elastic  moduli 
of  particulates  where  cementation  extends  far  beyond  immediate  grain-to-grain  contacts? 
Consider,  for  example,  a  close  random  pack  of  identical  glass  beads  where  the  pore  space 
is  completely  filled  with  epoxy.  This  aggregate  is  a  case  of  contacting  spherical  inclusions 
(glass)  embedded  in  the  epoxy  matrix.  From  CCT  we  know  that  the  filling  (epoxy)  located 
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at  the  contacts  contributes  most  to  the  stiffness  of  the  material.  Then  it  is  clear  that  a  direct 


use  of  an  inclusion  theory  (such  as  SC  or  DEM)  is  inappropriate  here  because  it  does  not 
account  for  the  special  properties  of  cemented  contacts.  Contact  cement  and  pore-filling 
cement  (Figure  Ic)  contribute  differently  to  the  stiffness  of  the  material. 


Figure  2.  CCT  theoretical  predictions  (curves)  and  experimental  data  (solid  circles)  for  the 
compressional  elastic  modulus  (the  product  of  the  bulk  density  and  P-wave  velocity  squared)  from 
elastic  wave  propagation  experiments.  The  cemented  samples  were  prepared  by  (a)  adding  epoxy  to  a 
close  random  pack  of  identical  elastic  spheres  (Yin,  1993),  and  (b)  adding  water  to  Ottawa  sand  and  then 
freezing  it  (Tutuncu  et  al.,  1997).  Both  epoxy  and  water  are  wetting  fluids  and,  when  they  solidify,  act 
as  contact  cement  in  the  form  of  pendular  rings.  The  modulus  is  plotted  versus  the  cement  saturation 
of  the  pore  space  of  the  uncemented  glass  bead  and  sand  packs. 

To  treat  this  problem  our  approach  is  to  obtain  a  solution  to  the  concentrated  inclusion 
problem  (the  bulk  and  shear  moduli  of  a  close  random  pack  of  identical  elastic  spheres 
whose  pore  space  is  partially  or  completely  filled  with  elastic  cement)  where  the  (previously 
ignored)  special  properties  of  the  contact-cement  part  of  the  matrix  are  recognized  and 
included.  The  problem  posed  is  important,  not  only  methodologically,  because  many 
cemented  geomaterials,  such  as,  e.g.,  asphalt  concrete,  have  inclusion  concentrations  close 
to  the  upper  limit,  small  void  concentration  (porosity)  and  large  amounts  of  pore-filling 
cement. 
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SOLUTION 


To  solve  this  problem,  we  assume  that  at  some  high  porosity  where  CCT  is  still 
applicable  (Figure  3,  point  B),  the  effective  bulk  and  shear  moduli  of  the  cemented  material 
are  those  of  a  homogeneous  isotropic  matrix  with  empty  spherical  inclusions  of  the  same 
porosity.  Next  we  use  an  effective  medium  theory,  for  example,  SC  or  DEM,  to  express 
these  effective  bulk  and  shear  moduli  through  the,  as  yet  unknown,  elastic  moduli  of  the 
hypothetical  matrix  material.  Then  the  elastic  moduli  of  the  matrix  material  can  be  found 
from  the  two  resulting  equations. 


Figure  3.  Proposed  solution  scheme.  "A",  uncemented  sphere  pack  of  zero  stiffness.  "B",  pack  with 
small  volume  of  contact  cement  whose  moduli  are  calculated  from  CCT;  the  moduli  of  the  hypothetical 
matrix  material  are  found  from  SC  or  DEM.  "C",  completely  cemented  matrix  and  sphere  pack.  "D", 
matrix  and  sphere  pack  of  non-zero  porosity.  The  graph  (compressional  modulus  versus  porosity)  is 
given  for  epoxy-cemented  glass  beads. 


To  find  the  elastic  moduli  of  the  sphere  pack  where  the  entire  pore  space  is  filled  with 
cement,  we  assume  that  they  are  the  same  as  those  of  the  hypothetical  matrix  whose 


inclusions  are  now  completely  filled  with  the  cement  (Figure  3,  point  C).  Those  moduli 
can  be  computed  using  the  chosen  effective  medium  theory. 

Finally,  we  find  the  moduli  of  the  cemented  sphere  pack  of  a  non-zero  porosity  by 
assuming  that  they  are  those  of  the  cement-filled  hypothetical  matrix  into  which  empty 
spherical  inclusions  are  embedded  (Figure  3,  point  D).  The  porosities  of  this  new  body 
with  inclusions  and  of  the  cemented  pack  are  the  same.  The  mathematical  expressions  for 
this  algorithm  are  given  in  Appendix. 

Of  course,  the  inclusions  do  not  have  to  be  spherical  (results  for  ellipsoidal-shaped 
inclusions  are  easily  obtained  within  the  same  approach),  and  effective  medium  theories 
other  than  SC  and  DEM  could  be  used. 


EXAMPLES 

In  the  examples  below,  the  compressional  (M),  bulk  (K),  and  shear  (G)  moduli  are 
calculated  from  the  measured  P-  and  S-wave  velocities  (V^  and  V^,  respectively)  and 
density  p  as  M  =  pV^,  G  =  pF/ ,  K  =  p{Vp  -  4F/  /  3). 

The  wave-length  in  all  the  experiments  is  much  larger  that  the  individual  grain  size, 
which  justifies  the  effective  medium  approach  to  the  problem.  Below,  we  use  the 
compressional  modulus  instead  of  the  bulk  modulus  because  the  former  is  calculated 
directly  from  data  which,  generally,  can  be  acquired  with  much  higher  accuracy  than 
data.  The  elastic  properties  of  materials  (glass,  quartz,  epoxy,  and  ice)  used  in  the 
calculations  are  given  in  Table  1. 


Table  1.  Elastic  moduli  of  materials  used  in  calculations. 


Material 

M  (GPa) 

G  (GPa) 

K  (GPa) 

Glass 

86 

30.34 

45.55 

Epoxy 

9.47 

2 

6.8 

Quartz 

96.6 

45 

36.6 

Ice 

13.27 

3.53 

8.57 
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Epoxy-Cemented  Glass  Beads.  The  elastic-wave  velocities  have  been  measured 
in  glass  bead  packs  cemented  with  epoxy  (Yin,  1993).  The  porosity  in  these  packs  varied 
from  zero  in  the  fully-cemented  pack  to  0.38  in  the  uncemented  pack.  In  its  liquid  state 
before  solidification  epoxy  is  a  wetting  fluid  and  tends  to  accumulate  at  the  grain  contacts. 
Therefore,  at  high  porosities  CCT  is  applicable  to  estimating  the  effective  elastic  moduli  of 
cemented  bead  packs.  We  use  this  theory  for  porosities  between  0.38  and  0.29.  For 
porosities  below  0.29  we  use  the  above-described  solution  with  SC  (Berryman,  1980)  and 
DEM  (Norris,  1985;  Zimmerman,  1991).  Inclusions  are  spherical. 

Our  theoretical  estimates  match  the  experimental  data  for  the  compressional  modulus 
extremely  well  (Figure  4)  and  overestimate  the  shear  modulus  (at  least  partially  because 
CCT  overestimates  the  shear  modulus  at  high  porosities).  We  cannot  judge  with  certainty 
whether  the  poor  prediction  of  the  shear  modulus  is  due  to  approximations  used  in  CCT  or 
to  systematic  errors  in  picking  the  5-wave  first  arrival  in  the  experiments. 


Figure  4.  Theoretical  (lines)  and  experimental  (solid  circles)  values  for  the  compressional  (a)  and  shear 
(b)  moduli  of  epoxy-cemented  glass  beads,  plotted  versus  porosity.  Bold  lines  are  from  CCT;  thin 
lines  are  from  the  new  solution. 
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An  adjustable  parameter  in  the  solution  is  porosity  above  which  CCT  is  used  and  below 
which  the  new  solution  is  employed.  The  new  solution  is  stable  with  respect  to  this 
parameter  (as  is  shown  in  Figure  5)  and,  therefore,  this  transition  porosity  can  be  arbitrarily 
chosen  within  the  reasonable  limits  where  CCT  is  applicable. 


Figure  5.  Stability  of  the  new  solution.  Compressional  and  shear  moduli  of  epoxy-cemented  glass 
beads  versus  porosity.  Bold  lines  are  from  CCT,  thin  lines  are  from  the  new  solution  (only  SC  is  used 
in  this  example).  They  start  at  a  varying  transitional  porosity.  Solid  circles  are  from  the  experiments. 

It  would  be  more  satisfying  if  all  of  these  curves  intersected  the  zero  porosity  axis  at  the 
same  modulus  value,  but  instead  we  find  that  the  present  formulation  of  the  theory 
produces  a  set  of  curves  that  are  shifted  (more  or  less)  in  parallel  depending  on  their  point 
of  intersection  with  the  CCT  curve.  This  result  suggests  that  the  current  formulation  might 
be  improved  using  a  more  sophisticated  effective  medium  method.  Further  discussion  of 
this  point  is  in  the  Conclusion  section  of  the  paper. 

Ice-Cemented  Ottawa  Sand.  The  elastic-wave  velocities  have  been  measured  in 
partially-saturated  and  frozen  pure-quartz  Ottawa  sand  (Tutuncu  et  al.,  1997).  Similar  to 
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the  first  example,  our  theoretical  solution  accurately  models  the  experimental  data  for  the 
compressional  modulus  and  overestimates  (as  CCT  does)  the  shear  modulus  (Figure  6). 


Figure  6.  Theoretical  (lines)  and  experimental  (solid  circles)  values  for  the  compressional  (a)  and  shear 
(b)  moduli  of  ice-cemented  Ottawa  sand,  plotted  versus  porosity.  Bold  lines  are  from  CCT,  thin  lines 
are  from  the  new  solution. 

The  high  degree  of  scatter  in  the  quartz-ice  data  (and  especially  so  for  the  shear 
modulus)  suggests  that  the  role  of  crystalline  ice  in  forming  the  microstructure  of  the 
cementing  contacts  between  quartz  grains  is  more  complex  (possibly  due  to  the  formation 
of  polycrystalline  ice)  than  in  the  other  examples  of  cemented  materials  considered  here. 

Sintered  Glass  Beads.  Elastic-wave  velocity  measurements  were  conducted  on 
sintered  glass  beads  by  Berge  et  al.  (1993).  Berryman  and  Berge  (1996)  discuss  the 
applicability  of  various  effective  medium  theories  to  this  dataset.  Microphotographs  show 
that  the  sintering  resulted  (because  of  glass  melting)  in  glass  pendular  rings  accumulated  at 
the  bead  contacts.  At  the  same  time,  microcracks  are  apparent  in  some  beads.  These 
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microcracks  are  likely  to  be  the  reason  for  the  data  scatter  observed  at  high  porosities 
(Figure  7).  Still,  the  high-porosity  data  can  be  reasonably  well  mimicked  by  CCT. 

The  SC-based  new  solution  matches  the  experimental  data  for  the  compressional 
modulus  extremely  well  (Figure  7a).  The  DEM  curve  (Figure  4)  underestimates  it.  Neither 
SC  nor  DEM  curves  match  the  shear  modulus  data  but  they  do  braeket  it.  In  this  case  we 
have  no  doubts  in  the  quality  of  the  data  because  the  datapoints  approach  the  (independently 
measured)  pure  glass  modulus  as  porosity  approaches  zero.  We  conclude  that  the  new 
solution  needs  to  be  improved  to  better  predict  the  shear  modulus  of  a  cemented  aggregate. 


Figure  7.  Theoretical  (lines)  and  experimental  (solid  circles)  values  for  the  compressional  (a)  and  shear 
(b)  moduli  of  sintered  glass  beads,  plotted  versus  porosity.  Bold  lines  are  from  CCT,  thin  lines  are 
from  the  new  solution. 

Natural  Quartz  Rocks.  Many  quartz  rocks  are  created  by  diagenesis  from  beach 
sand.  During  diagenesis,  quartz  rims  grow  on  the  surfaces  of  the  grains  cementing  them  in 
the  process.  This  cementation  scheme  is  different  from  the  pendular  ring  cementation 
because  part  of  the  cement  is  deposited  on  the  grain  surface  away  from  the  contacts.  At  the 
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same  porosity,  the  contact  cement  volume  and  thus  the  rock's  moduli  are  smaller  than  in  the 
pendular  ring  case.  Dvorkin  and  Nur  (1996)  apply  CCT  to  predicting  elastic  moduli  in 
high-porosity  sandstones.  The  results  match  the  data  (Strandenes,  1991)  well  (Figure  8). 

In  order  to  model  the  medium-  and  low-porosity  sandstone  data  (Han,  1986),  we  use 
CCT  with  the  pendular-ring  cementation.  The  datapoints  tend  to  approach  the  SC-based 
new  solution  line  as  porosity  approaches  zero  (Figure  8).  What  is  most  important,  the 
theoretical  results  at  zero  porosity  are  close  to  the  independently  determined  moduli  of  pure 
quartz.  The  DEM-based  new  solution  underestimates  these  moduli. 


Figure  8.  Theoretical  (lines)  and  experimental  (solid  circles)  values  for  the  compressional  (a)  and  shear 
(b)  moduli  of  natural  quartz  rocks,  plotted  versus  porosity.  Bold  lines  are  from  CCT,  thin  lines  are 
from  the  new  solution.  The  CCT  lines  are  plotted  for  two  cementation  schemes  (as  shown  in  the 
graphs):  the  pendular  rings  scheme  and  the  cementing  rims  scheme. 


CONCLUSION 

We  calculate  the  effective  moduli  of  an  elastic  matrix  with  elastic  inclusions  at  high 
concentration  by  modeling  it  as  a  cemented  particulate  aggregate.  This  approach  allows  us 
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to  account  for  the  important  role  of  the  contact  cement  part  of  the  matrix  and  address  the 
problem  of  directly  contacting  inclusions.  The  latter  problem  cannot  be  approached 
using  traditional  solutions  where  the  inclusions  are  considered  absolutely  rigid  (e.g., 
Frankel  and  Acrivos,  1967)  and  separated  by  a  finite-thickness  layer  of  the  matrix. 

The  solution  is  a  combination  of  the  contaet  eementation  theory  and  an  effective 
medium  theory  sueh  as  the  self-consistent  or  differential  effective  medium  approximation. 
The  solution  is  approximate  beeause  it  uses  a  number  of  assumptions.  Nevertheless,  its 
predictions  are  confirmed  by  experimental  data  on  cemented  glass  beads,  sand,  and  natural 
rocks.  The  self-consistent  approximation  gives  best  results  in  the  examples  considered. 
We  are  able  to  accurately  model  the  dynamic  compressional  moduli  of  these  aggregates.  At 
the  same  time,  the  errors  in  modeling  the  dynamic  shear  moduli  may  be  as  high  as  25%. 
There  may  be  at  least  two  reasons  for  these  errors:  (a)  commonly  encountered  inaeeuracy 
in  measuring  elastie  shear-wave  velocity;  and  (b)  theoretical  errors  in  modeling  the  effective 
elastic  modulus  of  a  cemented  high-porosity  sphere  pack.  For  the  latter  reason,  additional 
effort  is  needed  to  improve  our  theoretical  solution. 

One  legitimate  critieism  of  the  present  formulation  of  the  theory  is  that  neither  of  the 
two  effective  medium  theories  used,  although  well-established  and  relatively  easy  to  use, 
permits  the  user  to  control  the  microstructure  of  the  resulting  eomposite  in  any  substantial 
way.  In  particular,  both  approaches  when  applied  in  Steps  3  and  4  (see  Appendix)  will 
surely  violate  our  assumed  microstructure  for  intermediate  porosities  by  effectively  placing 
voids  in  the  grains  (where  we  do  not  want  them)  as  well  as  in  the  cement  material  (where 
we  really  do  want  them  to  be). 

This  effeet  will  be  comparatively  small  when  the  eementing  material  has  similar  or  the 
same  elastic  parameters  as  the  grains  as  is  the  case  for  our  glass-glass  and  quartz-quartz 
examples.  But  it  may  have  a  substantial  effeet  on  the  predicted  properties  in  the  eases 
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where  the  ratio  of  grain  to  cement  properties  is  large  as  is  the  case  for  glass-epoxy  and 
quartz-ice.  The  expected  bias  that  would  result  from  the  undesired  voids  in  the  grains  is 
that  the  predicted  effective  values  of  the  grain/cement  mixture  at  zero  porosity  must  contain 
a  higher  proportion  of  grain  material  in  the  total  volume  than  that  present  in  the  CCT 
calculations.  The  result  will  be  a  bias  towards  higher  than  realistic  values  at  zero  porosity 
and  this  is  exactly  what  we  observe  for  the  SC  theory  in  the  glass-epoxy  and  quartz-ice 
examples.  From  this  point  of  view,  the  present  results  should  be  considered  preliminary; 
we  need  to  continue  searching  for  relatively  simple  effective  medium  theories  that  will 
permit  us  to  gain  control  over  the  mierostructure  being  modeled. 

Still,  this  solution  can  be  directly  applied  to  estimating  the  stiffness  of  particulates  such 
as  asphalt  concrete  where  the  grains  touch  each  other  and  the  void  fraction  is  small. 
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APPENDIX:  SOLUTION  ALGORITHM 


Step  1,  Contact  Cementation  Theory  (CCT).  The  effective  compressional 
( ^ccr)-'  t’ulk  ( Kccr),  and  shear  ( Gca-)  moduli  of  a  random  dense  pack  of  identical  spheres 
cemented  at  their  contacts  are  (Dvorkin  and  Nur,  1996): 


V  _  0o)  e  ri  —^V  I  ^o)  f-<  o 

^CCT  —  r  '^CCT  —  r^ccr''  7^  ^c^f> 

6  5  20 


c’^T’  ^CCT  ~  ^CCT  2  ^CCT’ 


where  and  are  the  compressional  and  shear  moduli  of  the  cement  respectively; 
^0  ~  0.36  is  the  porosity  of  the  uncemented  pack;  and  n~S.5  is  the  average  number  of 
contacts  per  grain. 

Parameters  S„  and  are: 

=  A„(A„)a^  +  5„(A„)a  +  C„(AJ,  4(A„)  = -0.024153  •  A„-‘■^^^ 

5„(A„)  =  0.20405  •  C„(A„)  =  0.00024649  •  A„-' 

=  A,(A„  v)a^  +  5,(A„  v)a  +  C,(A^,  v), 

A,(A„  V)  =  -10"^  •  (2.26  +  2.07  v  +  2.3)  • 

V)  =  (0.0573v^  +  0.0937  v  + 0.202)  • 

C,(A„v)  =  10“^  -(9.654  +4.945v  +  3.1)-A,°°'*^’'''^®-^°" 

2a(i-v)a-v,) 

"  ttG  l-2v,  ’  nG 

where  G  and  v  are  the  shear  modulus  and  the  Poisson's  ratio  of  the  grain  material, 
respectively;  and  is  the  Poisson's  ratio  of  the  cement.  Parameter  a  depends  on  the 
contact  cementation  scheme.  It  is 

a  =  2[ 

3n(l-<Ao) 

for  the  pendular  ring  scheme  and 

^  ^  .2(^0 -0)  0,5 

3(1 -0o) 
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for  the  scheme  where  the  surface  of  a  grain  is  evenly  covered  with  the  cement  rim.  ^  is  the 
porosity  of  the  cemented  aggregate  which  is  smaller  than  because  of  cement  deposited  in 
the  pore  space  of  the  sphere  pack. 


Step  2,  Self-Consistent  Approximation  (SC).  If  the  moduli  and  of  the 

hypothetical  matrix  material  were  known,  then  the  self-consistent  approximation  for 
spherical  inclusions  of  void  with  concentration  ^  predicts  that  Kca-  and  G^cr  are  given  by 
the  coupled  equations  (Berryman,  1980): 
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1-^ 


^CCT  ^  ^CCT  ^  ^CCT  ^  ^CCT 
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1-^  ^  <t>  2  =  ^ccri^^ccr  +  ^G^cr) 


^CCT  Z  G^+Z  Z  +2G(’(^) 

These  equations  have  to  be  solved  iteratively  to  find  K^a-  and  G^cr  when  and  G^ 
are  known.  But,  since  and  G^cr  are  already  given  in  Step  1,  we  may  solve  the  two 
equations  instead  for  and  G^ .  Solving  in  this  direction,  the  equations  decouple  and 
provide  explicit,  unique  results  for  and  G^. 


Step  2',  Differential  Effective  Medium  Approximation  (DEM).  In  this  case 
and  G^  can  be  found  from  the  following  two  equations  (Norris,  1985): 

^ccr  _  /I  ^  ^s^CCT  ^  13  ^CCT  _  ^CCT  _ 0.75  -I-  _ 
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Norris'  result  can  be  rewritten  to  show  explicitly  that 
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where 


T  _  3(1  ^Vccr)  w  _  ^  ^^CCT  '^^CCT 
4{1+Vccr)'  2  3K,„  +  Gccr' 

Although  this  is  an  implicit  equation  relating  A,  to  (and  therefore  must  be  solved 
iteratively),  we  can  see  that  the  equation  is  "almost"  explicit  by  noting  that  “  0.2  is 
often  a  good  approximation  for  porous  media.  Then  both  A^  and  A^ct  will  be  small  in 
magnitude  and,  therefore,  the  last  factor  in  the  right-hand  side  of  the  equation  will  be 
closely  approximated  by  unity.  Thus  A^  is  uniquely  determined  by  Acer  0  in  the 
expected  range  of  the  elastic  parameters  for  porous  materials. 

Step  3  and  4,  SC.  The  moduli  Kp-n  and  Gp^  of  the  matrix  where  all  inclusions  (of 
concentration  0)  are  filled  with  the  cement  can  be  found  from  the  following  two  equations 
(Berryman,  1980): 

(1  -  0)(^,  -  Kp„)P^  +  0(^,  -  Kp,„)P,  =  0, 

(1  -  0)(G,  -  GpM  +  m  -  G,a,)Qc  =  0, 
where 

p  _  ^ Fill  ^  Gpjii  ^  _  ^Fill  “t*  ^  Gpiii 

K,  +  %G„'  K,+%G„' 

_  Gpiii  +  Z  ^  _  Gpii,  +  Z  ^  _  Gpjii(9Kpjii  +  8Gf,7/) 

G,  +  Z’  G,-hZ’  6(^,,,+2G„.„) 

Next  we  calculate  the  effective  moduli  ( Kp^  and  Gpff)  of  the  same  matrix  where  some 
of  the  inclusions  are  filled  with  the  cement,  and  some  are  empty.  If  the  concentration  of  the 
empty  inclusions  is  the  concentration  of  the  cemented  inclusions  is  because  the 

volumetric  fraction  of  the  matrix  is  still  1  -  0 .  The  desired  moduli  are  found  from  the 
following  two  equations: 

(1  -  -  Kpff)P,  +  -  Ke^Pc  -  tVo  =  0^ 

(1  -  <^)(G,  -  GpXQs  +  (0  -  -  Gp^Qc  -  (PeGsffQo  =  0, 
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where 
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Steps  3’  and  4%  DEM.  We  use  the  iterative  scheme  developed  by  Norris  (1985).  We 
choose  a  volumetric  fraction  increment  d(l>  such  that  d<j)«(f>.  The  total  number  of 
iterations  N  is  (Zimmerman,  1991): 
ln(l-0) 


N  =  -- 


d(j) 


The  number  of  the  first  steps  is  /  (j) .  During  these  iterations  we  calculate 


the  effective  bulk  and  shear  G^^  moduli  as: 
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The  last  N-N^  iterations  are  performed  using  equations 
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